Amplification of typhoon-generated near-inertial internal waves observed near the Tsushima oceanic front in the Sea of Japan

It is not fully understood how near-inertial kinetic energy (NIKE) is spatially distributed near Tsushima oceanic front (TOF) as a typhoon travels across the region. Underneath TOF, a year-round mooring covering a major part of water column was implemented in 2019. During summer, three massive typhoons (Krosa, Tapah, and Mitag) consecutively traversed the frontal area and delivered a substantial amount of NIKE into surface mixed layer. According to a mixed-layer slab model, NIKE was widely distributed near the cyclone’s track. The mooring observation exhibited the vertical distribution and pathways of surface-generated NIKE in response to the successive typhoon events. According to the modal decomposition, first three modes mostly explain the NIKE’s elevations following the typhoon events. According to ray-tracing experiments based on the internal-wave theory, large-scale near-inertial waves (NIWs) rapidly descend to a depth greater than 1000 m, while mesoscale NIWs slowly descend and rarely reached beyond the main pycnocline. Following the passage of Tapah, a profound energy mass was found nearly stationary at shallow depths coincident with vertical shear of geostrophic current. We infer that the descending rate of NIWs fell and then they were amplified through the energy conservation when the waves came from the north side of TOF.

www.nature.com/scientificreports/ captured near-inertial kinetic energy (NIKE) associated with a fast-travelling "bomb" cyclone over the sea during summer of 2015. Its NIKE was transferred to multiple inertial frequencies due to non-linear interactions amid the trapped waves. Intriguingly, the same cyclone in summer 2015 was also observed to cause profound turbulent mixing within another mesoscale eddy located in the opposite (i.e., Russian) side of the same ocean 16 . The propagation of NIKE can be restricted not only by the RVA body but also by vertical shear of the geostrophically balanced baroclinic current. That is, the lower propagation limit can be more strictly expressed as follows: where Ri g = N 2 ∂Ug ∂z 2 + ∂Vg ∂z 2 is the gradient Richardson number; N 2 = − g ρ ∂ρ ∂z is the buoyancy frequency; g the gravitational acceleration; z is vertical axis pointing upward; and ρ is the density of seawater, which represents the inclination of isopycnal interface 17 . According to their shipboard observations of internal waves and microscale turbulence in the TOF region, Kawaguchi et al. (2021) 18 demonstrated that the sloping isopycnals created geostrophic shear and limited the propagation of near-inertial waves toward the deeper layers. From a viewpoint of the wave kinetic energy distribution under the sea, it is worth evaluating the roles and impacts of the baroclinicity of the TOF as a function of the typhoon-induced near-inertial waves.
In the regions proximate to the Korean Peninsula, the previous studies revealed that the tidal forcing dominantly generated the internal waves and influenced the local turbulent mixing 19,20 . In the central part of Sea of Japan, the barotropic tidal currents are negligibly small 6,15,18,21 and are less likely subject to critical corruptions in the current data records by tides. In this sense, the TOF region is an ideal domain for the study of the interaction between oceanic front and wind-induced near-inertial waves.
Within the framework of the FRA-AORI Tsushima Warm Current Observatory (FATO) project, we implemented a year-round mooring observation of water current near the TOF in 2019-2020. To cover the entire depth where the surface generated near-inertial waves possibly travelled, we installed two sets of acoustic Doppler current profilers (ADCPs) and seven sets of single-point current meters (Fig. 1g) ("Data and methods"). In this paper, we analyzed the moored horizontal current with special focuses on the vertical distribution of NIKE, associated with the recurrent passages of massive typhoons. For the identification of the neighboring sources of near-inertial waves, we computed NIKE generated by the respective typhoon events near surface level across the entire Sea of Japan via the mixed layer slab model 3 . Based on the classical theory of the internal-wave dispersion 11 , the ray-tracing experiments were also conducted to interpret the distribution and vertical transfer of the observed www.nature.com/scientificreports/ signatures of NIKE. The satellite-based geostrophic currents were depicted as the supporting information of the TOF axis and neighboring mesoscale structures.

Results
Typhoon events in summer 2019. Between August and October, three massive typhoons, named Krosa, Tapah, and Mitag, successively or consecutively passed over the Sea of Japan and above the TOF (Table 1; Fig. 1d-f). The observed maximum wind speeds were 40, 35, and 40 m s −1 , while the minimum pressures near  the center were 965, 970, and 965 hPa, respectively. All the typhoons were associated with sufficiently warm surface temperatures in the tropical latitudes (10°-20°N) of the North Pacific Ocean. As for their trajectories, Krosa travelled through the middle of the Sea of Japan at a mean speed of 8.3 m s −1 . Tapah went nearly straightforward along the main island of Japan at 11.8 m s −1 . Mitag drew a crooked trajectory in clockwise direction through the southern part of the sea, at a mean travelling speed of 8.8 m s −1 .
Inertial wave kinetic energy input in surface mixed layer. The typhoons-delivered NIKE into SML was estimated from the slab model (see "Data and methods"). In Fig. 2, the cumulative NIKE input (Π) along with the cyclone's track is illustrated on the geographical map for each event.
The passage of Krosa produced a relatively small NIKE input, where the only noticeable patch was found on the right-hand side of the cyclone's track, centered at the geographical point (38°N, 136°E), whose magnitude was roughly Π = 7 kJ m −2 (Fig. 2a). Tapah delivered the maximum amount of NIKE at Π = 10-12 kJ m −2 , relatively evenly distributed at both sides of the cyclone's track (Fig. 2b). The maximum energy input was spotted roughly 100 km apart from the track, plausibly coincident to the core radius of the cyclone, where the tangential wind speed became greatest. Mitag produced an uneven distribution of the NIKE input (Fig. 2c). In the western half domain of the Sea of Japan, approximately along the cyclone track or slightly in its left-hand side, the greatest energy input can be found at the geographical point (39°N, 131-132°E), excessing Π = 15 kJ m −2 in the neighborhood. In the eastern domain, Mitag provided a relatively moderate NIKE input, in particular, at the right-hand side of its migration track 1 . The maximum surface energy input was found in the small region between the Noto Peninsula and the Sado Island, whose magnitude was approximately Π = 13-15 kJ m −2 .
At the nearest grid from the FATO station, the slab model simulations provided the accurate predictions of the inertial oscillation generation in the SML following the recurrent atmospheric events (Fig. 2d-h). The model predicted that there were outstanding oscillatory events in the SML in response to the three typhoon passages (Fig. 2e,f). During the passages of Krosa, Tapah, and Mitag, the SML was estimated to receive the net NIKE amounts of (i.e., dt ) of 5.5, 6.5, and 9 kJ m −2 , respectively (Fig. 2h). Following the three typhoon events, the slab model indicated that the SML inertial oscillation built up to be 0.5, 0.7, and 0.8 m s −1 , respectively, in the maximum amplitude of horizontal velocity (Fig. 2f). In addition to the typhoon events, the slab model showed that the transient gusty event that occurred in between Krosa and Tapah added a certain amount of NIKE (~ 2.0 kJ m −2 ) into the SML, generating an oscillation amplitude of 0.3 m s −1 (blue inverted triangle in Fig. 2).
Hydrographic and mesoscale features. Along the vertical profile of N(z) during T/S Oshoro-maru expedition, visiting the position of the FATO mooring, we distinguished the two distinct peaks in the stratification around depths of 30 and 200 m (Fig. 3a). The upper one indicated the base of SML and yielded N = 14 cycle per hour (cph), while the lower one indicated the main pycnocline and yielded N = 6 cph. At depths below the main pycnocline, N fell with depth until reaching nearly 0.2 cph right above the sea floor (~ 1800 m). The baroclinicity of the main pycnocline may play a crucial role in maintaining the geostrophic balance in the top layer above that pycnocline (Fig. 3).
According to the SSH spatial distribution, the main axis of near-surface current showed the prominent meandering feature near the FATO station (bold green in Fig. 1a-c). The maps depicted that the FATO site was located near the southern rim of a cyclonic mesoscale structure in August and September. Near the mooring station, the satellite image describes that the RVA was always positive for the whole summer months because of the cyclonic eddy, with an increase from ζ g = 0.05f to 0.2f (Fig. S1a).
For the longtime evolution of SSH and RVA horizontal distribution, we discover that the cyclonic eddy rapidly grew up in August and then remained intense throughout the following months of September and October ( Fig. S1a-g). Subsequently, the cyclonic feature stayed weak during the winter months, and eventually dissipated by February (Fig. S1h-j). www.nature.com/scientificreports/ Spectral signatures of near-inertial waves along vertical axis. The time-series observations of the horizontal current at the successive depths allowed for the spatial depiction of the vertical distribution of nearinertial oscillation intensity (Fig. 4). In the top layer (200 m deep or shallower), the sub-inertial frequencies of ω ≤ 10 -2 cph yielded a profound amount of kinetic energy, especially in the counterclockwise (CCW) rotational direction. This sub-inertial maximum was interpreted as a contribution from the surface-intensified geostrophic current, associated with the TOF jet. In (d-f), eastward and northward components are shown by red and blue, respectively. In (e) and (g), typhoon events of Krosa, Tapah, and Mitag are indicated by inverted triangles, respectively, in red, black, and white, and the gusty event in blue. www.nature.com/scientificreports/ According to the spectral analysis in the clockwise (CW) rotation, we could find a distinct peak around ω = f across the entire water column, including the bottom layer greater than 1000 m depth (Fig. 4a). In the top layer of the CW component, lying between the SML and the main pycnocline (hatched in Fig. 4), the near-inertial peak appeared to be relatively wide spreading at both sides of ω = f , i.e., 0.9-1.2f. When it went toward the deeper section, the peak appeared to be further sharpened, particularly at the sub-inertial side. The waves entering the deeper section were more strictly segregated based on the minimum propagation frequency of ω min = f due to no ambient flow leaving geostrophic shear significantly as well as the RVA terms in that depth range [Eq. (1)]. The near-inertial peak was found to be less significant in the counterclockwise (CCW) component than in the CW counterpart (Fig. 4b).
The energy contribution in the vicinity of ω = f (i.e., 0.9-1.2f) appeared to mostly result from the modulation locally near the FATO station, rather than from propagating waves from the neighboring latitudes, particularly from north (Fig. 4a). For example, the spectral energy can be found at ω = 1.1f , corresponding to the inertial frequency at 42.5°N, which is an unrealistic candidate as a source of typhoon-generated wave KE (Fig. 2a-c). www.nature.com/scientificreports/ The Doppler shifting ( �ω = k · U g ) could contribute to the wide-spreading frequencies around ω = f at the depths above the pycnocline (Fig. 4a). The large-scale waves associated with the cyclone speed can yield the simple estimate of the Doppler shift by �ω = |Ug | V cyc f . For example, providing the cyclone's travelling speed V cyc = 8 m s −1 (Table 1) and the magnitude of near-surface geostrophic current U g = 0.2-0.5 m s −1 would yield a frequency shift of 2-6% from the intrinsic wave frequency.
Wave properties deduced from oscillatory motions. From September 29 onward, the upward-looking ADCP at the top clearly demonstrated the periodic pattern of near-inertial horizontal velocity propagating in the upward direction over the depths of 100-350 m (Fig. 5a-c). According to the past studies (e.g., Leaman and Sanford 22 ), the upward propagation of phase accompanies oppositely directed (downward) propagation for wave kinetic energy. The best-fit regression line based on the phase evolution deduces a frequency of the mostdominant oscillation as ω = 1.05f , calculated at a depth of 153 m (Fig. 5d) (c.f., 1.02f in the equatorial Pacific 23 ).
Vertical shear of near-inertial current commonly depicts a clockwise rotation in hodographs (Fig. 5e). The trajectory was depicted with an hourly near-inertial band-passed velocity, for data chunks of 48 h with a halfoverlapping. We also drew best-fit ellipses for the hodographs by using the least-squares criterion technique 24 . When the near-inertial oscillation was the most intense, i.e., between September 29 and October 2, the rotational orbit was little compressed. The ratio between the semi-major and semi-minor axes, respectively a and b , was estimated as r = a/b = 1.1 . Aside from this, the ellipses were only slightly tilted from the horizontal by φ = 20° or smaller (Fig. 5f,g).
From the fact, the theory of WKB internal waves suggests that the packets of incident waves travelled approximately in the east-west direction 23 . Notwithstanding, we have only little confidence on the consequence because of the modest oblateness for the elliptical regressions (i.e., r ≈ 1 ). From the polarization viewpoint, the ratio r = a b can be related to that of the intrinsic frequency ( ω ) to the effective Coriolis frequency ( f e ) 23 : where k H = √ k 2 + l 2 and (k, l, m) are respectively zonal, meridional, and vertical wavenumbers, one can derive an estimate of horizontal wavelength ( H = 2π k H ) for the most dominant NIWs by an equation below: Here, we substituted the relation, C z = ω m , based on the vertical phase translation estimated from the ADCP observations (Fig. 5c). With measured ranges for the important parameters ( r = 1.1, ζ = 0.05-0.20 f 0 , N = 4-5 cph, and C z = 500-800 m d −1 ), we obtain the estimate of H falling within a typical value of 60-140 km. www.nature.com/scientificreports/ Horizontal wavelength and horizontal group velocity. The surface-delivered NIKE is expected to spatially radiate as internal waves. The horizontal migration is examined based on the rotational internal wave dispersion together with the principle vertical modes, n, as follows 25 : where k H is the horizontal wavenumber; c n is the eigenvalue phase speed of the n th vertical mode (i.e., n = 1, 2, 3, · · · ). The dynamic vertical mode calculation was performed using the averaged profile of N(z) during the T/S Oshoro-maru expedition in July ( The derivative of the wave dispersion (5) gives the horizontal group velocity ( Cg H(n) ) as the vertical mode solution for an arbitrary wavenumber: The analytic solution above illustrates the monotonous positive relationship between Cg H and k H , where it becomes constant to be Cg H = c n in the limit of k H = ∞ . It also suggests a linear increase of Cg H(n) with the wavenumber, for the bands of with R being a horizontal scale associated with the quasi-geostrophic balanced motions 12,26 . In the present case, R was estimated at roughly 100 km based on the mean profile of N(z) and its first baroclinic eigenvalue (Fig. 3). This is approximately consistent with the horizontal wavelength, 60-140 km, derived from the wave dispersion.
Kawaguchi and Yabe (unpublished) performed the high-resolving numerical simulations for the oceanic reactions to the wind forcing of Typhoon Tapah. They showed that the cyclone-resonance scale of 2π/k C H was the most profound irrespective of the presence or absence of the TOF. Also, they showed that the TOF resulted in the quasi-geostrophic disturbances, which empowered the near-inertial internal waves with the similar scale with them (i.e., R = 2π/k M H ). Regarding the mesoscale feature, the wave dispersion predicts the frequency at ω = 1.4f for the first baroclinic mode, which is obviously out of the typical range of the observed ω (Fig. 4). This implies that the first mode cannot exist as the propagating near-inertial waves, so that the second or greater vertical modes may instead determine the fastest travelling speed. Consequently, we estimated Cg H ≈ 20 km/day for both wavenumbers of k C H and k M H (Fig. 3e). It is necessary to clarify the potential influences of near-surface NIKE in remote places from the FATO station ( Fig. 2a-c). Overlain with the horizontal distribution of Π, the approximate radii were marked up over which near-inertial waves possibly travelled for a period of three and seven days (Fig. 2a-c). In general, the internal waves of k C H and k M H are commonly too sluggish to travel far horizontally and can reach ~ 60 km at most for three days. Therefore, it would not be significant at the FATO site immediately unless the typhoon supplies the inertial energy in its neighborhood. Here, we also maintain that the poleward propagation of near-inertial waves should be largely limited by the Coriolis term f (θ) in Eq. (1) as discussed in the past studies 27 .
Vertical distribution of wave kinetic energy. The combination of the observed horizontal current at many different depths clearly represented the vertical propagation of NIKE responding to the consecutive passages of the typhoons (Fig. 6). After the near-inertial bandpass filtering and the Wentzel-Kramers-Brillouin (WKB) normalization applied to the horizontal velocity (see the caption of Fig. 6; "Data and methods"), NIKE was dispersively transferred along the time-vertical dimension (Fig. 6c,d). In the late August and the early October, respectively, following Krosa and Tapah, a part of NIKE was apparently brought down to depths greater than 1000 m. Another part of NIKE appeared to remain in the vicinity of the lower pycnocline (e.g., 300 to 400 m). Following Tapah's passage in between the late September and the early October, there were extremely strong signals of near-inertial oscillations at depths of 100 m or shallower. The intensified NIKE remained at a similar depth for a few weeks and did not show any downward propagation until it dissipated eventually.
We also calculated the modal decomposition for near-inertial band-passed horizontal velocity (Fig. 6f). Please refer to "Data and methods" and Fig. S2 for more details. The decomposition clarifies the allocation of the observed NIKE by each vertical mode. In general, the first five vertical modes show the dominant contribution for the profound events in NIKE following typhoon passages, covering roughly 60-90% of the total amount. Regarding the events in August and early September, the first vertical mode explains roughly 35% of all, while the second and third modes do approximately 15 and 10%, respectively. During the greatest NIKE elevation observed following Tapah, the third mode claims the relatively large contribution by more than 40%, while the first and second modes share the remainder parts of 20-30%. www.nature.com/scientificreports/ To quantitatively discuss its vertical allocation, the wave-related NIKE was accumulated and compared between the two vertical segments of the top layer (50-300 m) and the bottom layer (300-1200 m), separated by the main pycnocline (Fig. 6e). For the estimate, we integrated the WKB-scaled NIKE over the N 0 -standardized vertical coordinate z * (see "Data and methods"). Overall, both cumulative curves displayed the similar responses to the cyclone passages. Commonly, the greatest peaks for the two segments lagged by about a week, relative to the surface energy input. The sum of NIKE was roughly comparable between the two segments. The peaks coincidental between the two layers suggested that a major part of the amount of wave kinetic energy was transported through the bottom by the near-inertial wave packet shortly after the cyclone's passage.
A notable exception of the near-inertial wave event was found to occur around October 1 (Fig. 6e). During this event, the maximum energy in the top layer, most likely as its response to Tapah, preceded the response of the bottom layer even by several days. Also, its integrated amount over the top layer reached 7 kJ m −2 or greater, which surpassed by far the response of the bottom layer (e.g., more than double). The maximum wave kinetic energy shortly attenuated to half or smaller in 10 days (Fig. 6e).
Following Tapah, the slab model estimated the gross amount of the surface NIKE input at around 6.5 kJ m −2 (Fig. 2h). Observed wave NIKE integrated over the top layer was comparable with the value estimated from the slab model. As the cause of excess NIKE of the top layer, we hypothesize the possibility of the temporary expansion of wave kinetic energy due to the rapid change in vertical group velocity. In the next section, we explore this hypothesis from the Lagrangian viewpoint for the near-inertial wave packet in the vertical direction.
As already mentioned above, NIKE in the upper part of water column showed a 5-6 days delay after the passage of Tapah (Fig. 6c-e) despite the slab model that simulates the immediate response to it in the SML (Fig. 2f). shear, (c) near-inertial band-passed (0.9-1.1f) current and (d) its magnitude. In (e), the vertical integral of NIKE is shown for depths of 50-300 m (red), 300-1200 m (blue) and the total (black). In (f), the contributions from the first five vertical modes are shown by shaded areas, with the total amount accumulating from n = 1 to 20 by bold blue curve, while the raw NIKE by bold black curve. Note that in (c-e) velocities are WKB-scaled (see "Data and methods" for more details). www.nature.com/scientificreports/ We hypothesize that a combination of potential factors caused the delayed response of NIWs and its noticeable amplitude. First, we focus on the sluggish geostrophic current near the surface that lasted for the period of September 21-28, which included the Tapah's passage (Fig. 6a). The sluggish current is likely due to the relative location of a cyclonic eddy (CE). We infer that the core of CE lied right over the station, resulting in such a weak geostrophic current observed. We suppose that the CE, characterized by a positive RVA, could block the NIW packets entering inside the structure, consistent with the numerous previous studies 5, 11,18,29 . During the shipboard survey near the FATO station in the late October, 2019, we indeed observed apparently weak signals of NIKE within the CE mesoscale structure compared to those in the ambient waters (Fig. 6 of Kawaguchi et al. 18 ). Second, we suggest a possible scenario that the delayed peak of NIKE was attributable for the late arrival of energy-containing NIW packets emanating from remote source regions. In fact, we spot the potential candidate as an abundant NIKE source, centered around (40°N, 136°E), at the distance of 5-day from the FATO station assuming the traveling speed of 20 km day −1 (Fig. 2b). For the third point, we know that vertical shear of the near-surface mean current concomitantly appeared along with the rapid growth of NIKE (Fig. 6b,c). From the fact, it is plausible that the mean current's vertical shear promoted the amplification of the NIW. This process will be discussed in the following section.

Discussion
We based our discussion on examining the paths of wave propagation, pursuing how it behaved depending on the characteristics of incident waves in response to the vertical group velocity ( Cg z ) as follows 11 : where the first term represents the vertical structure of density stratification N(z) ; the second and third terms are the cross product between vertical shear of geostrophic current  (Table 1) and at a depth of 50 m, with the initial conditions of (ω, k H ) and vertical shear of geostrophic ambient flow (U g , V g ) (Fig. 6b). The computation ended either at the maximum time of 60 d or when the particle arrived at a depth of 1000 m. In the ray-tracing experiments, the six different conditions of (ω, k H ) = 1.01f e , k C H , 1.03f e , k C H , 1.07f e , k C H , 1.01f e , k M H , 1.03f e , k M H , and 1.07f e , k M H were explicitly investigated (Fig. 7) (see "Data and methods" for further information).

Dependencies of wave rays on k H and ω.
In the default experiment, the shear term in Eq. (8) (see "Data and methods") was not included. Principally, the vertical distribution of N determines the common behaviors of the wave ray--the greatest descending speed can be viewed in the top layer, underneath the SML, followed by the significant attenuation near the main pycnocline at a depth of around 300 m (Fig. 7a). Then, the descending speed recovered in the deeper section underneath the pycnocline continuing through the bottom of the data record, 1200 m.
As was mentioned already, the horizontal wavenumber dominantly determined the magnitude of Cg z (Fig. 7a). The cyclone-internal wave resonant scale ( k H = k C H ) showed the greatest descending speed, typically travelling 1000 m within 5 d for any choice of ω . Meanwhile, the smaller scales of near-inertial waves ( k H = k M H ) generally showed sluggish downward energy propagation, resulting in the absolute postponement of the wave kinetic energy delivery to the deeper section. For the case of ω = 1.07f e , the mesoscale would take roughly 3 weeks or more. Accordingly, a variety of wave rays emanating from the consecutive typhoons of Tapah and Mitag intersected with each other, creating the intricate distribution of oscillatory wave kinetic energy across a major part of the water column.
In terms of the frequency of incident waves, it provided the modification in the descending speed of the wave packet (Fig. 7). Qualitatively, high ω tended to rapidly descend to get to the bottom, as opposed to low ω . Quantitatively, the wave packet of ω = 1.07f e arrived at 1200 m within 5-6 days, while ω = 1.01f e took at least 25 d. In fact, the observed distribution of NIKE was in close agreement with some of the predicted curves based on the first term of Eq. (8) (Fig. 7a). For example, following the passage of Krosa, the three distinctive wave rays were clearly visible. We can discover the guidance by the wave variants of ( ω , k H ) = (1.03f e , k M H ) descending to 900 m depth arriving around August 23, while (1.01f e , k C H ) and (1.03f e , k M H ) to 600 m depth around September 1 and 6, respectively.
As for Tapah, the variants of (1.01f e ,k C H ) and (1.03f e , k M H ) were only discernable, which may guide the energy down through the bottom of the current observation, 1200 m. For Mitag, the observed NIKE distribution was totally disturbed presumably due to the intersection of numerous rays from the consecutive typhoons of Tapah and Mitag. From the complicated NIKE distribution, any ray with significant amplitude was not discernible at the bottom layer beyond the main pycnocline.
Ray's refraction due to baroclinicity. Vertical shear of the baroclinic current provides a significant modification in the travelling path of the waves (Eq. (8)). Its intensity of the cross product ( ∂U g ∂z × k ) is maximized when an incident wavevector orthogonally meets the baroclinic flow (see the illustration in Fig. 8). Throughout summer and fall in 2019, the main TOF axis was oriented in the northeast or occasionally southeast directions (Fig. 1a-c). Vertical shear, ( ∂U g ∂z , ∂V g ∂z ), was horizontally in a similar direction with the surface geostrophic current since it is mostly surface intensified (Fig. 6a,b). The maximum magnitude of vertical shear was found around www.nature.com/scientificreports/ depths of 70-200 m. Baroclinic shear was generally noticeably profound in strength between the early August and the mid October as altering the primarily contributing components. The period of the maximum shear was fully overlapped with the recurrent typhoon passages (Table 1).
We examined the trajectories of the incoming waves emanating from the different source regions (south, west, north, and east) relative to the site (Fig. 7b-e). For Krosa, which delivered the greatest energy in the west and southwest (Fig. 2a), the mesoscale waves (e.g., 100 km) were significantly modulated to slow down as a result of vertical shear enhanced in the late September (Fig. 6c,d). Some rays from Krosa were persistent in the early October even after Tapah passed.
The rays from Tapah were relatively less modified than those from the others by the buoyancy term due to weakened vertical shear throughout October (Fig. 7). Exceptionally, the rays from north, in particular, with the mesoscale wavelengths showed the noticeable refraction beyond the buoyancy term. The downward-going www.nature.com/scientificreports/ waves appeared to bounce back toward the surface, coincident with the maximum of observed NIKE (Fig. 6e), and remained stagnant, overlapped with NIKE of the new-coming Mitag. Rays from Mitag were relatively less refracted by the shear term for any direction except for north (Fig. 2c). We maintain that near-inertial waves emanating from the south regions relative to the site is unlikely to be significant in amplitude due to the restriction of ω min (θ) for the northward propagation [i.e., Eq. (1)].

Concluding remarks
According to the mixed-layer slab model simulations, we depicted the horizontal distribution of surface inertial energy. During the passage of Krosa, the significant energy patch was positioned in the southwest direction about several days' distance from the mesoscale wavelengths (i.e., R ) (Fig. 2a). In our hypothesis, the southward component of the horizontal wave propagation could slow down the downward energy propagations. Tapah delivered the significant amount of NIKE within the accessible range, in particular, in the northern part of the FATO station. In particular, the mesoscale or smaller-scale wavelengths are more likely subject to the shear effect due to the baroclinicity, resulting in the modification of Cg z . In contrast, the cyclone-resonant scale is less likely subject to the effect because it immediately descends to the bottom layer before interacting with the front. Considering the conservation law of the wave action ( A = E ω , where E is the wave kinetic energy) with respect to the vertical direction, the rapid change of Cg z , resulting in the convergence, must be complemented by an increase of the wave amplitude 7 , i.e.: We here note that ω does not change along the trajectory 28 . This implies that the decelerated group velocity, associated with geostrophic shear, may result in the significant amplification of near-inertial waves. The phenomena based on the invariant action flux is also confirmed for near-inertial waves arrested within a negative RVA body 15,29,30 . In fact, near-inertial waves reacting with Tapah displayed the remarkable enhancement in the oscillation amplitude at a shallow depth of around 100 m (Fig. 6d). We suppose that the incident wave packet emanating from the northwestern energy source was enhanced through the amplification mechanism, as was mentioned above.
For the future perspectives, more thorough and more accurate estimate for the kinetic energy distribution associated with typhoon passages will require performing the in-situ observations in wider regions across the entire Sea of Japan.

Data and methods
The FATO mooring station was situated at the geographical point (38.72°N, 137.83°E), roughly 10 km northwest relative to the Sado Island (a mark in Fig. 1a-c) 15,18 . The full water depth at the site is 1770 m.

Current measurements.
In the present study, we mainly analyzed the horizontal current velocities for the period of July 1 and October 31, 2019. To track the vertical near-inertial wave propagation, we deployed a total of nine sets of current meters along the mooring system (Fig. 1g): 75 kHz and 300 kHz WH ADCPs (RD Instrument, US); three Infinity EMs (JFE Advantech Ltd., Japan); four 2 MHz AquaDopp 6000 m (Nortek, Norway) at respective nominal depths of 390, 420, 600, 700, 800, 900, 1000, 1100, and 1200 m. The 75 kHz and 300 kHz WH ADCPs were vertically oriented upward and downward, respectively. The top and bottom WH ADCPs measured (9) DA Dt = −A∇ · Cg ≈ −A ∂ Cg z ∂z Figure 8. Simplified schematic of ray refraction of near-inertial waves, associated with the vertical shear of TOF. The TOF current is assumed to be surface-intensified and go towards due east. When incident waves come from north as approaching the TOF, the ray cannot descend compared to the depth determined by the buoyancy term [Eq. (8)]. When waves come from south, they can go deeper than that. Please refer to "Data and methods" for further details. www.nature.com/scientificreports/ oceanic currents at 40 levels over a vertical range of 57-369 m and eight levels over 438-494 m, respectively, with a constant vertical interval of 8 m. Temporal interval was constantly 1 h for every current meter. The recorded current data were post-processed with a median filter for the removal of outliers. The estimated magnetic declination of 8.9 degrees in October 2019 (the World Magnetic Model) at the mooring location was corrected before the main analyses. The precision of horizontal velocity is 0.005 m s −1 for the RD and Nortek instruments and 0.01 m s −1 for the JFE Advantech instruments.
Hydrographic observations. During the shipboard observation of T/S Oshoro-Maru (Hokkaido University), we visited the FATO mooring station (more precisely at position within 1 nautical mile) on July 20 to take a couple of full-depth vertical profiles of conductivity-temperature-depth (CTD) variables by using an SBE-9plus (Seabird Electronics, Inc., US). The raw CTD profiles were averaged into a 1 dbar vertical pitch. The salinity was regularly calibrated with the high precision salinometer, AUTOSAL-8400B (Guildline Instruments, Canada) with the bottle sampled waters. The mean profile of N(z) is created as an average of the full-depth CTD profiling (Fig. 3a).
Surface geostrophic current from satellite altimetry. For the detection of mesoscale features near the FATO station, we used the satellite-based altimeter dataset distributed by the Copernicus Marine Environment Monitoring Service (CMES) (Fig. 1a-c). The sea surface height (SSH) and SSH-based geostrophic velocities were calculated during the multiple satellite investigations (e.g., TOPEX-Poseidon and Jason 1-3). The horizontal resolution was 1/4 degree (roughly 20 km). Temporal interval was daily. The geostrophically balanced RVA, ζ g , was computed at every grid cell, using the midpoint differential method 10 . Mixed layer slab model. For the source regions for the near-inertial wave packets by the typhoon events, the delivered wave kinetic energy input at the surface layer was calculated on the geographical map across the Sea of Japan. For this computation, we used the mixed layer slab model 3,31 . In the slab model, the kinetic energy of inertial motion in the SML is numerically solved thus: where Z I = U ML + iV ML is the inertial component of horizontal velocity; ρT = C D ρ a U a |U a | is the surface wind stress, where U a = U a + iV a is the 10-m height wind velocity in the complex form, and ρ a = 1.2 kg m −3 the air density. In the present study, we set the drag coefficient that is variable and dependent on the 10-m wind speed |U a | 32 : The wind-dependent coefficient yields the wind stress 29% and 57% greater for the cases of |U a | = 15 and 20 m s −1 if compared to those with the fixed value of C D = 1.3 × 10 −3 . We also add that the drag coefficient chosen here was widely used for the studies that quantify the impacts of atmospheric disturbances to the ocean 33,34 .
H is the mean thickness of SML, and we set it to be 20 m according to the hydrographic observation at the site (Fig. 3a). The superscript (*) signifies the conjugate of complex numbers, namely, D * = d − if and T * = τ x −iτ y ρ , where d is the damping constant for the inertial motion of the SML. We assumed in the calculation that d −1 = 4 day. The choice of d potentially brings some differences in the results. We evaluated its uncertainty by comparing the results for d −1 = 3-5 day and found a ~ 7% difference in the accumulated net NIKE flux at the surface (Fig. 2h) 35 .
For the calculation of the surface wind stress, the surface wind velocity was retrieved from the distributed dataset of "grid point value of mesoscale model (GPV-MSM)" published by the Japan Meteorological Agency (JMA) (Fig. 1d-f). The horizontal resolution was 1/20° in latitude and 1/16° in longitude. The temporal interval was hourly.
We stress that the amount of wave kinetic energy supplied by a travelling typhoon is not necessarily dependent of its intensity. Alternatively, its travelling speed can be another important factor to determine the oscillation strength. In other words, a slow-moving cyclone may be subject to the cancellation by the preexisting counterflow in the same layer. The fast-moving cyclone has a reduced chance of the cancellation.

Ray-tracing experiments.
For an understanding of the observed energy allocation, it provides a merit to track the pathways of near-inertial wave migration in the vertical sense. We conducted the ray-tracing investigation based on the simplified dispersion relationship of near-inertial waves 11 . Along the vertical axis, the location of a certain wave packet was tracked by integrating the vertical group velocity in time thus: In the continuously stratified water column, we assume 11 : www.nature.com/scientificreports/ where k = (k, l, m) is the wavenumber vector; and k H = √ k 2 + l 2 is the magnitude of the horizontal wavenumber components. For the simplicity, the vertical wavenumber was determined using the following simplified version of the wave dispersion 36 : For the range of ω = 1.01-1.07f, one would obtain the rough estimates of 1500-4000 m for the cycloneresonant scale of 680 km ( = 2π/k C H ), while 250-650 m for the mesoscale disturbances of 100 km ( = 2π/k M H ). In the calculation above, we used the canonical value of N = N 0 = 3 cph 37 . Dividing the full depth of ~ 1800 m by the estimated vertical wavelength would give a simple estimate of the vertical mode; k C H corresponds to n = 1 , while k M H does the higher modes of n ≥ 3. Equation (13) was re-arranged as follows: The above relationship describes the dependency of Cg z on the horizontal scale of propagating waves--with a large-scale internal wave, the descending rate of a wave packet tends to increase, as opposed to a small-scale internal wave. The frequency ω and horizontal wavenumber k H are explicitly given to examine the influences on the results.
In this study, we considered the frequencies of ω = 1.01f e , 1.03f e , and 1.07f e as the feasible range of the propagating waves. For the horizontal wavenumber, we tested the following options of k C H = 9 × 10 -6 rad m −1 and k M H = 6 × 10 -5 rad m −1 , respectively, representing the scales associated with the cyclone-induced resonance and the mesoscale features. For N(z) , we used the 1-m average from the onsite expedition in July (Fig. 3a). Vertical shear of geostrophic current, ( ∂U g ∂z , ∂V g ∂z ), was obtained as the vertical gradient of the 24-h low-passed horizontal velocity (see Fig. 6b). Subsequently, it was interpolated to the locations of N(z) data points along the vertical coordinate.
In Eq. (13), the first term mostly represents the vertical structure of the density stratification, referred to as the buoyancy term. Meanwhile, the second term represents vertical shear of the geostrophic current, referred to as the shear term, with a vertical component of the cross product with the horizontal wavenumber vector k , namely, We considered some realistic situations that an internal wave packet approached the geostrophic jet that flowed in parallel to the TOF (Fig. 8). On the horizontal plane, the directions of the phase and group velocities were identical. If a wave packet comes from the left side relative to the jet (i.e., ∂U ∂z l − ∂V ∂z k < 0 ) the downward migration of the packet slows down or is flipped in the principal direction toward the surface (Fig. 8). If the wave comes from the right ( ∂U ∂z l − ∂V ∂z k > 0 ), it boosts the downward speed of the wave packet migration.
Application of WKB scaling. For the fair representation of near-inertial wave amplitude, the near-inertial band-passed horizontal velocities U NI = (U NI , V NI ) were rescaled based on the WKB theory 24 by using the observed vertical profile of N(z) thus (Fig. 3a): where N 0 = 3 cph representing the canonical intensity of the density stratification 37 . For the integrating procedure of the WKB-scaled wave kinetic energy, the vertical coordinate was stretched by: The new vertical coordinate yields the maximum depth of about 750 m (c.f., 1200 m in the original coordinate).
Modal decomposition. The modal decomposition of the observed horizontal velocity is performed with the least-squared technique based on the eigenfunctions derived from the in-situ profile of N(z) 38 (Fig. 6f). For the vertical profile, we used the mean profile obtained during the Oshoro-maru cruise in July (Fig. 3a). The arrays of the observed near-inertial band-passed current (Fig. 6c) are linearly interpolated and extrapolated onto the uniformly spaced vertical coordinate at the 1-m interval, between the depths of 0 and 1200 m, corresponding to the original CTD profile.
The modal solutions describe the zonal motion of seawater for each vertical mode, n , at an arbitrary depth: www.nature.com/scientificreports/ Regarding the meridional velocity, v , it is defined and solved in the similar way with u . Here, we assume the WKB-type solutions (Fig. 3b,c)

expressed by:
where h is the bottom depth. The equation assures the boundary condition of φ n = 0 at z = 0 and z = −h . In the present study, we applied the boundary condition at the greatest observational depth of 1200 m. The average value of N is given by Using the normal modes and the near-inertial band-passed horizontal velocity from the mooring observation, ( U n , V n ) are estimated by solving the following matrix equations: where the vectors, A, B, and , are given by where j signifies the integer number of the vertical grids with J being its largest grid number; M is the greatest number of the vertical mode. In the equation above, we used the simplified representation for the derivative function: φ n = ∂φ n ∂z .

Data availability
Output data of the GPV reanalysis model developed by the Japan Meteorological Association can be downloaded from the website below: http:// datab ase. rish. kyoto-u. ac. jp/ arch/ jmada ta/ data/ gpv/ origi nal. The rest of dataset is available upon reasonable request. Correspondence and request for the materials should be addressed to Y.K.